from sklearn.preprocessing import MinMaxScaler
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, fcluster
import pandas as pd
import numpy as np


def AHCRBF(datafile, trainNum, resultfile):
    # 读数据划分训练集
    oriData = np.array(pd.read_excel(datafile))
    colNum = oriData.shape[1]-1   # 获取输入层结点数
    oriDataTrain = oriData[0:trainNum, :]
    data = oriDataTrain[:, 0:colNum]
    totalOutput = oriDataTrain[:, colNum].reshape(trainNum, 1)

    # 归一化
    scalerx = MinMaxScaler(feature_range=(-1, 1))
    xNormData = scalerx.fit_transform(data)
    scalary = MinMaxScaler(feature_range=(-1, 1))
    yNormData = scalary.fit_transform(totalOutput)

    # 初始化参数
    inputNums = colNum
    maxcount = 1000
    samplenum = data.shape[0]
    precision = 0.001
    alpha = 0.005  # 学习率设定值
    hideNums = 9
    error = np.zeros((1, maxcount+1))
    w = np.random.random((hideNums, 1))  # w表示隐层到输出层的权值

    # 获取距离矩阵，第二参数指定距离计算方法
    yy = pdist(xNormData, 'euclidean')
    # AHC_RBF_alogrithm.m:37
    # 获取聚类，第二参数指定层次聚类方式
    zz = linkage(yy, 'median')
    # 指定获取簇类个数
    Ncluster = hideNums
    # 获取指定Ncluster个数的聚类结果
    c = fcluster(zz, Ncluster, 'maxclust').reshape(trainNum, 1)
    AGNES_location_calculation_matrix = np.concatenate((xNormData, c), axis=1)

    xNormdata = xNormData
    Idx = np.zeros((samplenum, 1))
    Idx = c

    count = np.zeros((hideNums, 1))
    sum1 = np.zeros((hideNums, inputNums))
    for i in range(samplenum):
        if 1 == Idx[i]:
            count[0] = count[0] + 1
            for j in range(inputNums):
                sum1[0, j] = sum1[0, j] + xNormdata[i, j]
        elif 2 == Idx[i]:
            count[1] = count[1] + 1
            for j in range(inputNums):
                sum1[1, j] = sum1[1, j] + xNormdata[i, j]
        elif 3 == Idx[i]:
            count[2] = count[2] + 1
            for j in range(inputNums):
                sum1[2, j] = sum1[2, j] + xNormdata[i, j]
        elif 4 == Idx[i]:
            count[3] = count[3] + 1
            for j in range(inputNums):
                sum1[3, j] = sum1[3, j] + xNormdata[i, j]
        elif 5 == Idx[i]:
            count[4] = count[4] + 1
            for j in range(inputNums):
                sum1[4, j] = sum1[4, j] + xNormdata[i, j]
        elif 6 == Idx[i]:
            count[5] = count[5] + 1
            for j in range(inputNums):
                sum1[5, j] = sum1[5, j] + xNormdata[i, j]
        elif 7 == Idx[i]:
            count[6] = count[6] + 1
            for j in range(inputNums):
                sum1[6, j] = sum1[6, j] + xNormdata[i, j]
        elif 8 == Idx[i]:
            count[7] = count[7] + 1
            for j in range(inputNums):
                sum1[7, j] = sum1[7, j] + xNormdata[i, j]
        elif 9 == Idx[i]:
            count[8] = count[8] + 1
            for j in range(inputNums):
                sum1[8, j] = sum1[8, j] + xNormdata[i, j]

    # 计算聚类中心的坐标
    average = np.zeros((hideNums, inputNums))
    for i in range(0, hideNums):
        for j in range(0, inputNums):
            average[i, j] = sum1[i, j] / count[i]
    C = average

    # 求扩展常数
    dd = np.zeros((1, hideNums))
    for i in range(0, hideNums):
        dmin = 10000
        for j in range(0, hideNums):
            dist = 0
            for k in range(0, inputNums):
                dist = dist + (C[i, k] - C[j, k]) ** 2
            if dist < dmin and i != j:
                dmin = dist
        dd[0, i] = dmin

    # b为隐含层的输出向量Hj
    b = np.zeros((samplenum, hideNums))
    for i in range(0, samplenum):
        for j in range(0, hideNums):
            sum1 = 0
            for k in range(0, inputNums):
                sum1 = sum1 + (xNormData[i, k] - C[j, k]) ** 2
            b[i, j] = np.exp(- sum1 / (2*dd[0, j]))

    count = 0
    while count < maxcount:
        c = 0
        errorp = np.zeros((1, samplenum))
        while c < samplenum:
            o = 0.0
            for j in range(0, hideNums):
                o = o + b[c, j]*w[j, 0]
            # 反馈/修改;
            errortmp = (yNormData[c] - o) ** 2
            errorp[0, c] = errorp[0, c] + errortmp
            errorp[0, c] = 0.5*errorp[0, c]
            yitao = yNormData[c] - o
            for j in range(0, hideNums):
                w[j, 0] = w[j, 0] + alpha*yitao*b[c, j]
            c = c + 1
        # 求最后一次迭代的误差
        tmp = 0.0
        for i in range(0, samplenum):
            tmp = tmp + errorp[0, i]*errorp[0, i]
        tmp = tmp / c
        error[0, count] = np.sqrt(tmp)
        if error[0, count] < precision:
            break
        count = count + 1

    # ---------------------------------------------------
    # 预测
    rowNum = oriData.shape[0]
    oriDataTest = oriData[trainNum:rowNum, :]
    dataTest = oriDataTest[:, 0:colNum]
    sampleTest = dataTest.shape[0]
    theoryOutputTest = oriDataTest[:, colNum]
    yDataTest = theoryOutputTest
    xNormDataTest = scalerx.transform(dataTest)

    errRel = np.zeros((sampleTest, 1))
    avgErrRel = 0.0
    RMSE = 0.0
    absHideOut = np.zeros((sampleTest, hideNums))
    netActualOut = np.zeros((sampleTest, 1))
    for i in range(0, sampleTest):
        for j in range(0, hideNums):
            sum1 = 0
            for k in range(0, inputNums):
                sum1 = sum1 + (xNormDataTest[i, k] - C[j, k]) ** 2
            absHideOut[i, j] = np.exp(- sum1 / (2*dd[0, j]))
        absOut = 0.0
        for r in range(0, hideNums):
            absOut = absOut + absHideOut[i, r]*w[r, 0]
        netActualOut[i] = absOut

    netOut = scalary.inverse_transform(netActualOut)

    for i in range(0, sampleTest):
        if yDataTest[i] != 0:
            errRel[i] = errRel[i] + abs((yDataTest[i] - netOut[i])) / yDataTest[i]
        else:
            errRel[i] = errRel[i] + abs(netOut[i])
        RMSE = RMSE + (yDataTest[i] - netOut[i]) ** 2
        avgErrRel = avgErrRel + errRel[i]

    avgErrRel = np.mat(avgErrRel)*100 / sampleTest  # 矩阵形式
    avgErrRel = avgErrRel[0, 0]
    RMSE = RMSE / sampleTest
    RMSE = np.sqrt(RMSE)
    RMSE = RMSE[0]

    netOut = pd.DataFrame(netOut)
    netOut.to_excel(resultfile, index=None, header=None)

    return [avgErrRel, RMSE]